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Abstract. Following on after three previous papers discussing the formation of 
primordial black holes during the radiation-dominated era of the early universe, we 
present here results for a related toy problem where the parameter w in the equation 
of state is allowed to take values different from 1/3. This can give useful insight into 
the main case of interest. We consider also the possibility of having a cosmological 
constant, interpreted as a vacuum energy component (corresponding to w = — 1), 
revising the results presented in the first paper of this series with use of the present 
more-refined version of our numerical code. As in our previous work, we have started 
our simulations with initial supra-horizon scale perturbations of a type which could 
have come from inflation, characterized by a pure growing- mode component. We have 
verified that the critical collapse behaviour is preserved in this context for all of the 
different values of w considered, matching well with the parameter values found for 
the scaling laws in more idealized circumstances. 
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1. Introduction 

A population of black holes might have been formed in the early universe by means 
of gravitational collapse of cosmological perturbations corresponding to perturbations 
of the metric with large-enough amplitude (Zel'dovich & Novikov (1969) pQ; Hawking 
(1971) [2]; Carr [31 H]). Various groups have then demonstrated this numerically, finding 
that these primordial black holes (PBHs) would have been formed if some cosmological 
perturbations had reached an amplitude 5, at the time of cosmological horizon crossing, 
greater than a certain threshold value 5 C (Nadezhin, Novikov & Polnarev (1978) [5]; 
Bicknell & Henriksen (1979) [5J; Niemeyer and Jedamzik (1999) [7J; Shibata and Sasaki 
(1999) [8]; Hawke & Stewart (2002) [9.]; Musco, Miller & Rezzolla (2005) [10J). This 
5 is often defined as the relative mass excess inside the overdense region, measured at 
the time of horizon crossing. The simulations performed in these studies clarified many 
different aspects of the nature of the collapse with particular reference to the radiative 
era of the universe (with the equation of state being taken as p = |e, where p is the 
pressure, e is the energy density and c 2 s = 1/3 is the sound speed squared). This context 
represents the most commonly studied scenario for PBH formation. 

In 1999, Niemeyer & Jedamzik [TJ [TT] showed that the masses of PBHs produced in 
the radiative era by perturbations with a given profile type, follow the typical scaling- 
law behaviour of critical collapse, first discovered for idealized conditions by Choptuik 
[13], i.e. the masses of the black holes produced follow a power law Mbh oc (5 — <5 C ) 7 if 

5 is close enough to 8 C . More generally, one can get a critical collapse for a perfect fluid 
when the equation of state considered has the form p = we, with w varying between 
and 1, as shown by Neilsen & Choptuik in 2000 [14] using a succession of "imploding 
shells" of matter as the initial conditions. Moreover, they found that the value of the 
critical exponent 7 was dependent only on the value of w and not on the particular form 
of the perturbation profile. 

In 2002, Hawke & Stewart [9] came back to the problem of critical collapse in the 
context of PBH formation in the early universe with a radiative fluid, and investigated 
the nature of the collapse going down to smaller values of (5 — 5 C ) than Niemeyer 

6 Jedamzik had been able to do with their code. For the larger values of (5 — S c ), 
comparable to those of [7], Hawke & Stewart again found a scaling law with a similar 
value of 7, but for smaller (5 — 5 C ) they saw formation of strong shocks and the curve 
flattened off at a minimum mass of around 10~ 3 of their horizon mass. They also found 
that the value of S c depended strongly on the shape of the perturbation, in contrast to 

In our 2005 paper [10], we considered PBH formation with initial conditions 
given by linear perturbations of the energy density and velocity fields approximating 
a cosmological growing-mode, and imposed with a length-scale much larger than the 
cosmological horizon. The domination of the growing component becomes even greater 
by the time of horizon crossing, as any residual decaying component dies away. We again 
saw a scaling law with similar 7 for the range of values of (5 — 5 C ) used by [7] but found 
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very different values of 5 C from [7] for similar perturbation shapes (we were concentrating 
on those types of perturbation). This difference was attributed to the fact that in [7J, 
the perturbations were made just in the energy density and were imposed directly at 
the horizon scale, so that their value of S c was calculated with inclusion of a substantial 
decaying component which did not then contribute to the black hole formation. If one 
focuses on the effect of perturbations originating from inflation, then it is clearly only 
growing-mode- type perturbations which are relevant at much later times. 

In our 2007 paper [15], this kind of cosmologically-relevant initial perturbation was 
imposed in a more precise way by using an asymptotic quasi-homogeneous solution 
[16J. Starting from a curvature perturbation, which is a time- independent quantity 
when the perturbation length-scale is much larger than the cosmological horizon [17] . 
perturbations in all of the other quantities can then be specified consistently, giving a 
solution which behaves as a pure growing mode in the linear regime. It is possible in 
this way to discriminate very clearly between cosmological perturbations that will give 
rise to PBHs and ones that will just disperse into the surrounding medium. We found 
that only when there is a large-amplitude perturbation of the metric, corresponding to 
a non-linear initial perturbation of the curvature, can one get 5 > 5 C at horizon crossing, 
giving rise to a gravitational collapse leading to black-hole formation. 

In our 2009 paper [TH], we returned to the issue raised by Hawke & Stewart [H] 
concerning whether the scaling law would continue down to very small values of (5 — 5 c ). 
In order to address this, we needed to modify our code with the inclusion of adaptive 
mesh refinement (AMR) so as to be able to handle the extreme conditions which arise 
near to the critical limit. In view of our earlier work, we were particularly wanting to 
investigate the effect of using growing-mode-type initial perturbations imposed outside 
the horizon scale, rather than the non-linear sub-horizon scale initial perturbations used 
in [§]. We again used the quasi- homogeneous solution to provide our initial conditions. 
Doing this, we did not observe the shock formation seen in [9J (which we attributed to 
the presence of a non-linear decaying component in their calculations) but instead got a 
regular scaling law behaviour all the way down to the vicinity of the resolution limit of 
our scheme (going beyond the most extreme values shown in [9]). A striking feature of 
our calculations was the appearance of an "intermediate state" during collapse for cases 
near to the critical limit: the collapse proceeded to a compactness 2M/R ~ 0.5 (where 
R and M are the current radius and mass of the condensation) which then remained 
roughly constant for a large number of dynamical timescales, with the condensation 
progressively shedding matter via a relativistic wind and shrinking, until eventually it 
either reached a mass at which the final collapse to a black hole could take place or 
started to disperse. 

The earlier work focused on taking w = 1/3, which is the most interesting case 
from a cosmological point of view. However, for gaining more insight, it is useful to 
consider also other values for w, in order to get a broader view. Here we consider a 
range of different values, following exactly the same procedures for problem setup and 
calculation as used previously for w = 1/3. We want to see whether the scaling law 



Critical collapse and black hole formation within an expanding perfect fluid. 



4 



behaviour continues to be present, as in the ideal cases [Hj, and to analyze how the 
values of 5 C and 7 change with w, as well as clarifying the nature of the intermediate 
state. Another aim is to revisit our previous calculation in [TO] regarding critical collapse 
within a radiative fluid in the presence of a cosmological constant (considered as a 
vacuum energy giving a fluid component with w = —1). The results presented in this 
paper has been obtained using the same numerical code as in [18] , with some fine tuning. 

After this introduction, Sections 2 and 3 give brief summaries of our mathematical 
formulation of the problem, and of the numerical methods being used for the simulations; 
Section 4 contains the results which we have obtained from investigating the effects of 
varying w and from revisiting the case including a cosmological constant; Section 5 
contains conclusions. Throughout, we use units for which c = G = 1, except where 
otherwise stated. 

2. Mathematical formulation of the problem 

For the calculations described here, we have followed the same basic methodology as 
described in Papers 1, 2 and 3 (respectively [TOJ [151 HH])- We therefore give just a brief 
summary of it here; more details are contained in the previous papers. 

We use two different formulations of the general relativistic hydrodynamic 
equations: one for setting the initial conditions and the other for studying the black 
hole formation. Throughout, we are assuming spherical symmetry and that the medium 
can be treated as a perfect fluid; we use a Lagrangian formulation of the equations with 
a radial coordinate r which is co-moving with the matter. 

For setting the initial conditions, it is convenient to use a diagonal form of the 
metric, with the time coordinate t reducing to the standard FRW time in the case of a 
homogeneous medium with no perturbations. (This sort of time coordinate is therefore 
often referred to as "cosmic time"). We write this metric in the form given by Misner 
Sz Sharp (MS) [19J, whose approach we follow here in writing the GR hydrodynamic 
equations: 



with the coefficients a, b and R being functions of r and t and with R playing the role 
of an Eulerian radial coordinate. Using the definitions 



rfs 2 = _ a 2 df 2 + b 2 dr 2 + R 2 ( M 2 + ^2 ^2) 



(1) 




(2) 




(3) 



one defines the quantities 
U = D t R, 



(4) 



and 



r = 



D r R, 



(5) 
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where U is the radial component of four-velocity in the "Eulerian" frame and T is a 
generalized Lorentz factor. The metric coefficient b can then be written as 

b =f^ (6) 

In general p represents the rest mass density but, when the fluid is relativistic, as 
for a radiative medium, the rest-mass makes a negligible contribution to the energy 
density and p represents just the compression factor. With this specification the GR 
hydrodynamic equations can then be written in the following form (where we use the 
notation that e is the energy density, p is the pressure and M is the mass contained 
inside radius R): 



DM 



M 

D r p + — + A-nRp 



L(e+p) ' R 2 



(7) 



D t p = -^D r (R 2 U), (8) 

D t e = m-Dtp, (9) 

D t M = -4irR 2 pU, (10) 
a 

D r a = D r p, (11) 

e + p 

D r M = 47iR 2 Te, (12) 

plus a constraint equation 

F 2 = i + u 2 - — . (13) 
R 

An equation of state is also needed, and we are here considering ones of the form 

p = we (14) 

where if is a constant parameter which can take various values. When a cosmological 
constant is added, the total energy density e and the total pressure p are given by 

A A 
e = e+— p = p -— (15) 

where A/87T is the vacuum component of the energy density and —A/8n is the vacuum 



component of the pressure, which are related by equation (14) with w — — 1 



The initial conditions are set by introducing a perturbation of the otherwise uniform 
medium representing the cosmological background solution, with the length-scale of the 
perturbation R being much larger than the cosmological horizon Rh = H -1 . Under 
these circumstances, the perturbations in e and U can be extremely small while still 
giving a large-amplitude perturbation of the metric (as is necessary if a black hole 
is eventually to be formed) and the above system of equations can then be solved 
analytically to first order in the small parameter e = (Rh/Rq) 2 « 1. A full discussion 
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of this has been given in Paper 2 and the result, referred to as the "quasi-homogeneous 
solution" , gives formulae for the perturbations of all of the metric and hydrodynamical 
quantities in terms only of a curvature perturbation profile K(r) that is conveniently 
time- independent when e << 1 |17j . 

To characterize the amplitude of the perturbation, we use the integrated quantity 

which measures the relative mass excess within the overdense region, as frequently done 
in the literature; e& is the background value of the energy density and ro is the co-moving 
length-scale of the overdense region of the perturbation, determined by the particular 
expression used for K(r). In the present work, as already done in Papers 2 and 3 (to 
which we refer for details), we make the simplest choice of using a Gaussian for the 
initial curvature profile K(r), which leads to a corresponding mexican hat profile for 
the energy density at the initial time: 

/ q /„ \2\ 

(17) 



cxp 





e b (t) 3(l + «/) 

(19) 

where A = ro/v^B is the independent parameter used at initial time to vary the 
perturbation amplitude. It can then be shown that the time variation of 5 is given 
by 

5 + 3w A 2 ( 3 



^ = ^W^f^VV (20) 

within the linear regime where e << 1. It follows that 5 then has the familiar growth 
with time||j 

2(l + 3m) 
/ I \ 3(l+») 

<t) = <U) [jj • (21) 

The MS approach using cosmic time slicing is convenient for setting initial 
conditions but has the well-known drawback for calculations of black hole formation 
that singularities are formed rather quickly when using it and then further, potentially 
observable, evolution cannot be followed unless an excision procedure is used. Various 
slicing conditions can be used to avoid this difficulty but for calculations in spherical 
symmetry it is particularly convenient to use null slicing. In our work we have used 

| To be strictly precise, as explained in detail in Papers 2 and 3, we evaluate 8 analytically at the 
initial time taking into account only the time-independent part of the energy density perturbation 
profile, corresponding to e = 1. This gives a value of 6 of the same order as the one that would be 
measured at horizon crossing time. 
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the "observer time" null-slicing formulation of Hernandez & Misner [20J where the time 
coordinate is taken as the time at which an outgoing radial light ray emanating from an 
event reaches a distant observer. (In the original formulation, this observer was placed 
at future null infinity but for calculations in an expanding cosmological background we 
use an FRW fundamental observer sufficiently far from the perturbed region so as to 
be unaffected by the perturbation.) We use the MS approach for setting up initial data 
and for evolving it so as to produce data on a null slice. This is then used as input for 
our observer-time code with which we follow the black hole formation. 
For the observer-time calculation, the metric ([I]) is re- written as 

ds 2 = -f 2 du 2 - 2fb dr du + R 2 (d6 2 + sin 2 Od<p 2 ) , (22) 

where u is the observer time and / is the new lapse function. The operators equivalent 
to g and Q are now 

Dt = U^-), (23) 





( d_ 


7 








5 


V dr 



D ^A^)' (24) 

where D k is the radial derivative in the null slice and the corresponding derivative in 
the Misner-Sharp space-like slice is given by 

D r = D k - D t . (25) 

The hydrodynamic equations can then be formulated in an analogous way to what was 
done in cosmic time. We will not repeat this here but refer the reader to our discussion 

in nm. 



3. The method of calculation 

The present calculations have been made with the same code as that used for Paper 3. 
Since this has been fully described previously, we will just give a brief outline of it here. 
It is an explicit Lagrangian hydrodynamics code based on that of Miller & Motta (1989) 
[2T] but with the grid organized in a way similar to that of Miller & Rezzolla (1995) 
[22] which was designed for calculations in an expanding cosmological background. The 
code has a long history and has been carefully tested in its various forms. The basic 
grid uses logarithmic spacing in a mass-type comoving coordinate, allowing it to reach 
out to very large radii while giving finer resolution at small radii. Our initial data is 
derived from the quasi-homogeneous solution and is specified on a space-like slice (at 
constant cosmic time) with e = 1CT 2 , giving Rq = 10 R H . The outer edge of the grid 
has usually been placed at 90 Rh, which is far enough away so that there is no causal 
contact between it and the perturbed region during the time of the calculations. The 
initial data is then evolved using the MS equations ((7 13), so as to generate a second set 



of initial data on a null slice (at constant observer time) and the null-slice initial data 
is then evolved using the observer-time equations [TP"] . 
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For the calculations presented in Paper 3, we introduced an adaptive mesh 
refinement scheme (AMR), on top of the existing logarithmic grid, giving us sufficient 
resolution so as to be able to follow black hole formation down to extremely small values 
of (5 — 5 C ). Having the AMR is particularly important for allowing us to follow the deep 
voids which form outside the central collapsing region in cases very close to the critical 
limit. The same scheme has been used with just minor modifications for the calculations 
presented here. Our aim in writing the AMR was to avoid the use of artificial viscosity, 
but it has now emerged that some residual artificial viscosity was still present and is 
necessary for correct functioning of the code in its present form. The presence of this 
is not thought to affect the results presented in any important way, however. We have 
successfully used the scheme with more than thirty levels of refinement and all relevant 
features of the solutions have been fully resolved. 



4. Description of the results 

Here we describe results from numerical simulations of black hole formation in an 
expanding medium, generalizing those obtained for PBH formation in the radiative 
era of the early universe by using values for the constant w in the equation of state 
different from 1/3. This is to be considered as a toy problem, but one which can give 
useful insight into the real cases of interest. We are presenting results for w in the 
range between and 0.6. As mentioned above, nearly critical collapse to form black 
holes within this type of expanding background is characterized by a relativistic wind 
forming outside the central condensation which excavates a deep void around it. The 
wind and the void become more extreme for a given (5 — 5 C ) as w is increased, becoming 
very challenging for our AMR scheme, even though it has gone to more than thirty 
levels of refinement. We therefore limit ourselves to a maximum w of 0.6. 

As we know from the theory of cosmological perturbations, a growing supra-horizon- 
scale overdensity in the early universe increases in both amplitude and length-scale until 
it enters within the cosmological horizon, at which point it can start to collapse if the 
density amplitude reached is large enough to counteract the cosmological expansion. 
However, the internal pressure of the fluid works against the collapse, tending to disperse 
the condensation into the surrounding medium. Whether or not the collapse successfully 
proceeds to completion is determined by whether the perturbation amplitude 5 is large 
enough, and the value of 5 is directly connected with the initial curvature perturbation 



as shown by expressions (17), (18) and (20). The competition between gravity and the 



pressure gradients determines the threshold amplitude 5 C that discriminates between 
perturbations which are able to form PBHs (5 > 5 C ) and perturbations that disperse 
into the surrounding medium (5 < 5 C ). An approximate estimate for the value of this 
parameter, which is crucial for determining the abundance of PBHs that could be formed 
in the early universe, was made for the first time in 1975 in the landmark paper of Carr 
[3], where from considerations based on the Jeans length it was found that 

5 C ~ w . (26) 
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Figure 1. Dependence of 5 C on the value of w. The dotted line shows Carr's 
estimate for 5 C ~ w based on a Jeans-length argument. 



However a precise evaluation can be obtained only with the aid of numerical simulations, 
solving the non-linear system of general relativistic hydrodynamical equations. We have 
focused on the case of a radiative fluid with w = 1/3, finding 5 C ~ 0.45, but in general 
5 C would vary if w is allowed to vary, as Carr already noticed, because it measures the 
relative contributions of energy density and pressure in the fluid. 

In figure [TJ we present results obtained with our code for the dependence of 5 C on 
w, showing also the comparison of the numerical results with Carr's estimate, 5 C ~ w, 
which is plotted as a dashed line. One can see that w — 1/2 represents a dividing 
line: 5 C > w for w < 1/2 while 5 C < w for w > 1/2. This relation between 5 C and w 
confirms that any epochs in the early universe when the equation of state softens would 
be favourable for enhancing PBH production. One example of this could be the QCD 
phase transition, which has previously been discussed in connection with a PBH model 
for MACHOs |12j . and there could also be some possibility of PBHs formed towards 
the end of the radiative era (when their masses could be higher) providing seeds for 
supermassive black holes. However any realistic discussion of this would need to take 
into account many other physical processes in addition to just a softening of the equation 
of state. 

As previously mentioned, nearly-critical collapse is characterized by the appearance 
of a long-lived intermediate state for the forming condensation, where a close balance 
between gravity and the mass-shedding caused by the surface pressure gradients 
maintains an almost constant value of its compactness (as given by the peak value of 
2M/R) which can last for many dynamical timescales. The same behaviour is observed 
for each of the different values of w, with the value of (2M / 'R) pea k varying with w in a 
similar way to <5 C ; indeed, these two quantities have approximately the same value (that 
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Figure 2. The behaviour of 2M/R is plotted against R/Rh for a nearly- 
critical case (5 ~ S c ) for w = 0.1 and w = 1/3, at different time levels. Rh is 
the cosmological horizon scale at the moment of horizon crossing. The dashed 
curve shows the initial conditions used by the observer-time code. 



of S c being slightly smaller). In figure |2j we show two examples of the intermediate 
state, for w — 0.1 and w = 1/3 (the other cases follow the same pattern), plotting the 
behaviour of 2M/R against R/Rh at different time levels (Rh being the cosmological 
horizon scale at the moment of horizon crossing). The dashed curve shows the initial 
conditions used by the observer-time code, after which the evolution proceeds in the 
direction of decreasing (2M / 'R) pe ak until the intermediate state is reached. There is 
then a roughly self-similar behaviour, with the profile of 2M/R moving inwards with 
an almost constant peak value. Because increasing the value of w increases the ratio 
between the internal pressure and the energy density, it is consistent that higher values 
of w correspond to a higher value of [2M / R) pea k for the intermediate state. 

In figure [3] we show the scaling law behaviour of Mbh with respect to (5 — S c ) for 
several values of w ranging from 0.1 to 0.6. Good scaling laws are obtained in each case 
over the range of values of (5 — 5 C ) that we were able to study (note that this range was 
much smaller in the cases with w > 1/3 because of the extreme conditions encountered 
there). The best-fit values of 7 are perfectly consistent with those obtained by Neilsen 
and Choptuik [T3] in their imploding shells calculation. We have also observed that 7 
itself scales roughly linearly with w, as one can see from figure [4], where we have used 
the results for 7 obtained in [T3] for the values of w > 0.6 which we were not able to 
calculate ourselves. The dashed line represents the straight-line best fit to the data. 

In Paper 1 we investigated critical collapse for a radiative fluid in the presence of 



a cosmological constant A, as described in section 2 here by equations (14) and (15). 



Results from that calculation are shown again here as the left-hand panel of Figure [5] 
and it is noticeable that the scaling laws with and without inclusion of a finite A have 
significantly different gradients 7. At the time of Paper 1, the code was not able to get 
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Figure 3. Scaling behaviour for Mbh as function of (5 — 5 C ) for values of w 
ranging from 0.1 to 0.6. Mu is the cosmological horizon scale at the moment 
of horizon crossing. The best-fit values for 7 are indicated in each case. 



very close to the critical limit because of not yet having the AMR, and the lowest value 
of (5 — 5 C ) that we were able to treat was around 10~ 3 . We have now made a similar 
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Figure 4. The behaviour of 7 as a function of w, showing the rather good 
approximation given by the straight-line fit 7 = 0.81u> + 0.09. 

calculation with the present version of the code, going down to much smaller values of 
(5 — S c ), and the results are shown in the right-hand panel of Figure [5] On this scale, 
the difference between the gradients seen previously is hardly visible (the region of the 
left-hand panel is shown by the dashed lines in the upper right-hand corner of this plot). 
The difference is only present for values of 5 not very close to S c , with Mbh > O.IMh, 
while for smaller values of Mbh the densities within the condensation are higher and 
the contribution from the non-zero A becomes negligible by comparison, so that the 
value for 7 becomes just that for the fluid alone. We can then identify two regimes 
for the scaling law with A > 0: one where (A/87re) is mainly much less than 1 during 
the collapse, corresponding to the lower part of the scaling law with the same value 
of 7 as for A = 0, and the other where (A/87re) ~ 1 during the collapse, giving the 
changed value for 7. These two regimes can be identified respectively in the upper-right 
and lower-left parts of the right-hand panel of Figure |5j separated approximately by 
the dashed line. One can notice also a small offset between the masses even when the 
effect of A during the collapse is negligible. This is due to the dependence on A of the 
cosmological horizon mass at the moment of horizon crossing, Mjj, which we are using 
as our unit of mass. 

5. Conclusions 

Following on after our previous work investigating primordial black hole formation 
during the radiative era of the early universe, we have here presented results for a 
related toy problem where the equation of state parameter w is allowed to take constant 
values different from the radiative value of 1/3, the aim being to gain further insight into 
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Figure 5. Scaling behaviour for Mbh as function of (6 — 8 C ) for a radiative 
perfect fluid (w = 1/3), with and without a cosmological constant. The left- 
hand panel reproduces Figure 1 from Paper 1 |10| . while the right-hand panel 
shows results from the new calculation with the same value of A, made with 
the current version of the code. The region of the left-hand panel is shown by 
the dashed lines in the upper right-hand corner of the right-hand panel. 



our main case of interest by considering it within a broader context. Our calculations 
have been particularly focused on examining how the critical collapse phenomenon 
is manifested for primordial black hole formation in the early universe, initiated by 
perturbations of a type which could have come from inflation, and we follow self- 
consistently both the collapse to form the black hole and the continuing expansion 
of the universe. We use a purpose-built Lagrangian AMR code which allows a very 
large dynamic range. 

In this paper, we have followed formation of black holes for values of w ranging 
from 0.1 to 0.6 and, for the type of perturbation that we are studying, we find that 
scaling-law behaviour persists down to the smallest masses that we are able to follow. 
The critical threshold amplitude 5 C , the intermediate state compactness (2M/R) peak and 
the scaling-law exponent 7 all vary with w in a way which is consistent with the fact 
that w measures the ratio between the pressure and energy density of the fluid. In 
particular, the values obtained for the exponent 7 are in almost perfect agreement with 
the results obtained by Neilsen and Choptuik [13] in their imploding shells calculation. 
We have also observed that 7 and w follow a roughly linear relationship. 

In addition we have revisited the results of Paper 1 concerning scaling-law behaviour 
for a radiative fluid in the presence of a cosmological constant. We have seen that the 
there are two regimes of the scaling law, with different values of 7, corresponding to 
whether or not A/87T remains comparable with the fluid energy density during the 
collapse. 
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